Exercise for Recurrent Neural Networks: PdM Regression


By Prof. Seungchul Lee
http://iai.postech.ac.kr/
Industrial AI Lab at POSTECH

1. Problem Definition¶

  • Predictive Maintenance
    • This predictive maintenance project focuses on the techniques used to predict when an in-service machine will fail, so that maintenance can be planned in advance.


  • Overview of the RUL estimation strategy
    • The remaining life of a test unit is estimated based on the actual life of a training unit that has the most similar degradation pattern.


  • Data description
    • Training data: It is the aircraft engine run-to-failure data.
    • Testing data: It is the aircraft engine operating data without failure events recorded.
    • Ground truth data: It contains the information of true remaining cycles for each engine in the testing data.



  • Operational conditions
    • Altitude (0 ~ 42,000 feet), Mach number (0 ~ 0.84) and throttle resolver angle (20 ~ 100)
  • 21 sensors
    • Different measurements related to the engine state at runtime.

2. Import Library¶

In [1]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
from sklearn import preprocessing

3. Load Data¶

In [2]:
train_df = pd.read_csv('./data_files/PM_train.txt', sep=" ", header=None)
train_df.head()
Out[2]:
0 1 2 3 4 5 6 7 8 9 ... 18 19 20 21 22 23 24 25 26 27
0 1 1 -0.0007 -0.0004 100.0 518.67 641.82 1589.70 1400.60 14.62 ... 8138.62 8.4195 0.03 392 2388 100.0 39.06 23.4190 NaN NaN
1 1 2 0.0019 -0.0003 100.0 518.67 642.15 1591.82 1403.14 14.62 ... 8131.49 8.4318 0.03 392 2388 100.0 39.00 23.4236 NaN NaN
2 1 3 -0.0043 0.0003 100.0 518.67 642.35 1587.99 1404.20 14.62 ... 8133.23 8.4178 0.03 390 2388 100.0 38.95 23.3442 NaN NaN
3 1 4 0.0007 0.0000 100.0 518.67 642.35 1582.79 1401.87 14.62 ... 8133.83 8.3682 0.03 392 2388 100.0 38.88 23.3739 NaN NaN
4 1 5 -0.0019 -0.0002 100.0 518.67 642.37 1582.85 1406.22 14.62 ... 8133.80 8.4294 0.03 393 2388 100.0 38.90 23.4044 NaN NaN

5 rows × 28 columns

In [3]:
test_df = pd.read_csv('./data_files/PM_test.txt', sep=" ", header=None)
test_df.head()
Out[3]:
0 1 2 3 4 5 6 7 8 9 ... 18 19 20 21 22 23 24 25 26 27
0 1 1 0.0023 0.0003 100.0 518.67 643.02 1585.29 1398.21 14.62 ... 8125.55 8.4052 0.03 392 2388 100.0 38.86 23.3735 NaN NaN
1 1 2 -0.0027 -0.0003 100.0 518.67 641.71 1588.45 1395.42 14.62 ... 8139.62 8.3803 0.03 393 2388 100.0 39.02 23.3916 NaN NaN
2 1 3 0.0003 0.0001 100.0 518.67 642.46 1586.94 1401.34 14.62 ... 8130.10 8.4441 0.03 393 2388 100.0 39.08 23.4166 NaN NaN
3 1 4 0.0042 0.0000 100.0 518.67 642.44 1584.12 1406.42 14.62 ... 8132.90 8.3917 0.03 391 2388 100.0 39.00 23.3737 NaN NaN
4 1 5 0.0014 0.0000 100.0 518.67 642.51 1587.19 1401.92 14.62 ... 8129.54 8.4031 0.03 390 2388 100.0 38.99 23.4130 NaN NaN

5 rows × 28 columns

In [4]:
truth_df = pd.read_csv('./data_files/PM_truth.txt', sep=" ", header=None)
truth_df.head()
Out[4]:
0 1
0 112 NaN
1 98 NaN
2 69 NaN
3 82 NaN
4 91 NaN

3.1. Plot¶

In [5]:
def plot_sensor_data(data, purpose):
    for j in data.keys().difference([0,1,26,27]).values:

        max_cycle = data.groupby([0]).count()[j].max()
        
        if purpose == 'train':
            plt.figure(figsize = (20, 1.5))
            for i in range(data[0].max()):
                padding = max_cycle - data[data[0] == i][j].values.shape[0]
                plt.plot(np.arange(padding, max_cycle), data[data[0] == i][j].values)
            if j == 2 or j == 3 or j == 4:
                plt.title(str(j-1) + 'th setting')
            else:
                plt.title(str(j-4) + 'th sensor')
            plt.vlines(max_cycle, data[j].min(), data[j].max(), colors = 'r', linestyles = 'dashed')
            plt.text(max_cycle+1, data[j].median(), "failure", fontsize=14)
            plt.show()
            
        if purpose == 'test':
            plt.figure(figsize = (20, 1.5))
            for i in range(data[0].max()):
                plt.plot(data[data[0] == i][j].values)
            if j == 2 or j == 3 or j == 4:
                plt.title(str(j-1) + 'th setting')
            else:
                plt.title(str(j-4) + 'th sensor')
            plt.text(max_cycle+1, data[j].median(), "on working", fontsize=14)
            plt.show()
In [6]:
# plot train data
plot_sensor_data(train_df, 'train')
In [7]:
# plot test data
plot_sensor_data(test_df, 'test')

4. Pre-processing¶

4.1. Drop NaN column¶

In [8]:
train_df.drop(train_df.columns[[26, 27]], axis=1, inplace=True)
train_df.columns = ['id', 'cycle', 'setting1', 'setting2', 'setting3', 's1', 's2', 's3',
                     's4', 's5', 's6', 's7', 's8', 's9', 's10', 's11', 's12', 's13', 's14',
                     's15', 's16', 's17', 's18', 's19', 's20', 's21']


test_df.drop(test_df.columns[[26, 27]], axis=1, inplace=True)
test_df.columns = ['id', 'cycle', 'setting1', 'setting2', 'setting3', 's1', 's2', 's3',
                     's4', 's5', 's6', 's7', 's8', 's9', 's10', 's11', 's12', 's13', 's14',
                     's15', 's16', 's17', 's18', 's19', 's20', 's21']

truth_df.drop(truth_df.columns[[1]], axis=1, inplace=True)

4.2. Train Data¶

In [9]:
train_df.head()
Out[9]:
id cycle setting1 setting2 setting3 s1 s2 s3 s4 s5 ... s12 s13 s14 s15 s16 s17 s18 s19 s20 s21
0 1 1 -0.0007 -0.0004 100.0 518.67 641.82 1589.70 1400.60 14.62 ... 521.66 2388.02 8138.62 8.4195 0.03 392 2388 100.0 39.06 23.4190
1 1 2 0.0019 -0.0003 100.0 518.67 642.15 1591.82 1403.14 14.62 ... 522.28 2388.07 8131.49 8.4318 0.03 392 2388 100.0 39.00 23.4236
2 1 3 -0.0043 0.0003 100.0 518.67 642.35 1587.99 1404.20 14.62 ... 522.42 2388.03 8133.23 8.4178 0.03 390 2388 100.0 38.95 23.3442
3 1 4 0.0007 0.0000 100.0 518.67 642.35 1582.79 1401.87 14.62 ... 522.86 2388.08 8133.83 8.3682 0.03 392 2388 100.0 38.88 23.3739
4 1 5 -0.0019 -0.0002 100.0 518.67 642.37 1582.85 1406.22 14.62 ... 522.19 2388.04 8133.80 8.4294 0.03 393 2388 100.0 38.90 23.4044

5 rows × 26 columns

4.2.1. Labeling¶

In [10]:
train_df['RUL'] = train_df.groupby(['id'])['cycle'].transform(max) - train_df['cycle']
train_df.head()
Out[10]:
id cycle setting1 setting2 setting3 s1 s2 s3 s4 s5 ... s13 s14 s15 s16 s17 s18 s19 s20 s21 RUL
0 1 1 -0.0007 -0.0004 100.0 518.67 641.82 1589.70 1400.60 14.62 ... 2388.02 8138.62 8.4195 0.03 392 2388 100.0 39.06 23.4190 191
1 1 2 0.0019 -0.0003 100.0 518.67 642.15 1591.82 1403.14 14.62 ... 2388.07 8131.49 8.4318 0.03 392 2388 100.0 39.00 23.4236 190
2 1 3 -0.0043 0.0003 100.0 518.67 642.35 1587.99 1404.20 14.62 ... 2388.03 8133.23 8.4178 0.03 390 2388 100.0 38.95 23.3442 189
3 1 4 0.0007 0.0000 100.0 518.67 642.35 1582.79 1401.87 14.62 ... 2388.08 8133.83 8.3682 0.03 392 2388 100.0 38.88 23.3739 188
4 1 5 -0.0019 -0.0002 100.0 518.67 642.37 1582.85 1406.22 14.62 ... 2388.04 8133.80 8.4294 0.03 393 2388 100.0 38.90 23.4044 187

5 rows × 27 columns

4.2.2. Normalize¶

In [11]:
cols_normalize = train_df.columns.difference(['id','cycle','RUL'])
cols_normalize
Out[11]:
Index(['s1', 's10', 's11', 's12', 's13', 's14', 's15', 's16', 's17', 's18',
       's19', 's2', 's20', 's21', 's3', 's4', 's5', 's6', 's7', 's8', 's9',
       'setting1', 'setting2', 'setting3'],
      dtype='object')
In [12]:
# normalize
min_max_scaler = preprocessing.MinMaxScaler()
norm_train_df = pd.DataFrame(min_max_scaler.fit_transform(train_df[cols_normalize]),
                             columns=cols_normalize,
                             index=train_df.index)

norm_train_df.head()
Out[12]:
s1 s10 s11 s12 s13 s14 s15 s16 s17 s18 ... s3 s4 s5 s6 s7 s8 s9 setting1 setting2 setting3
0 0.0 0.0 0.369048 0.633262 0.205882 0.199608 0.363986 0.0 0.333333 0.0 ... 0.406802 0.309757 0.0 1.0 0.726248 0.242424 0.109755 0.459770 0.166667 0.0
1 0.0 0.0 0.380952 0.765458 0.279412 0.162813 0.411312 0.0 0.333333 0.0 ... 0.453019 0.352633 0.0 1.0 0.628019 0.212121 0.100242 0.609195 0.250000 0.0
2 0.0 0.0 0.250000 0.795309 0.220588 0.171793 0.357445 0.0 0.166667 0.0 ... 0.369523 0.370527 0.0 1.0 0.710145 0.272727 0.140043 0.252874 0.750000 0.0
3 0.0 0.0 0.166667 0.889126 0.294118 0.174889 0.166603 0.0 0.333333 0.0 ... 0.256159 0.331195 0.0 1.0 0.740741 0.318182 0.124518 0.540230 0.500000 0.0
4 0.0 0.0 0.255952 0.746269 0.235294 0.174734 0.402078 0.0 0.416667 0.0 ... 0.257467 0.404625 0.0 1.0 0.668277 0.242424 0.149960 0.390805 0.333333 0.0

5 rows × 24 columns

In [13]:
join_df = train_df[train_df.columns.difference(cols_normalize)].join(norm_train_df)

train_df = join_df.reindex(columns = train_df.columns)

train_df.head()
Out[13]:
id cycle setting1 setting2 setting3 s1 s2 s3 s4 s5 ... s13 s14 s15 s16 s17 s18 s19 s20 s21 RUL
0 1 1 0.459770 0.166667 0.0 0.0 0.183735 0.406802 0.309757 0.0 ... 0.205882 0.199608 0.363986 0.0 0.333333 0.0 0.0 0.713178 0.724662 191
1 1 2 0.609195 0.250000 0.0 0.0 0.283133 0.453019 0.352633 0.0 ... 0.279412 0.162813 0.411312 0.0 0.333333 0.0 0.0 0.666667 0.731014 190
2 1 3 0.252874 0.750000 0.0 0.0 0.343373 0.369523 0.370527 0.0 ... 0.220588 0.171793 0.357445 0.0 0.166667 0.0 0.0 0.627907 0.621375 189
3 1 4 0.540230 0.500000 0.0 0.0 0.343373 0.256159 0.331195 0.0 ... 0.294118 0.174889 0.166603 0.0 0.333333 0.0 0.0 0.573643 0.662386 188
4 1 5 0.390805 0.333333 0.0 0.0 0.349398 0.257467 0.404625 0.0 ... 0.235294 0.174734 0.402078 0.0 0.416667 0.0 0.0 0.589147 0.704502 187

5 rows × 27 columns

4.3. Test data¶

4.3.1. Labeling¶

In [14]:
# Get value of max cycle
rul = pd.DataFrame(test_df.groupby('id')['cycle'].max()).reset_index()
rul.columns = ['id', 'max']
rul.head()
Out[14]:
id max
0 1 31
1 2 49
2 3 126
3 4 106
4 5 98
In [15]:
truth_df.columns = ['more']
truth_df['id'] = truth_df.index + 1
truth_df.head()
Out[15]:
more id
0 112 1
1 98 2
2 69 3
3 82 4
4 91 5
In [16]:
truth_df['max'] = truth_df['more'] + rul['max']
truth_df.drop('more', axis=1, inplace=True)
truth_df.head()
Out[16]:
id max
0 1 143
1 2 147
2 3 195
3 4 188
4 5 189
In [17]:
# Get RUL
test_df = test_df.merge(truth_df)

test_df['RUL'] = test_df['max'] - test_df['cycle']
test_df.drop('max', axis=1, inplace=True)

test_df.head()
Out[17]:
id cycle setting1 setting2 setting3 s1 s2 s3 s4 s5 ... s13 s14 s15 s16 s17 s18 s19 s20 s21 RUL
0 1 1 0.0023 0.0003 100.0 518.67 643.02 1585.29 1398.21 14.62 ... 2388.03 8125.55 8.4052 0.03 392 2388 100.0 38.86 23.3735 142
1 1 2 -0.0027 -0.0003 100.0 518.67 641.71 1588.45 1395.42 14.62 ... 2388.06 8139.62 8.3803 0.03 393 2388 100.0 39.02 23.3916 141
2 1 3 0.0003 0.0001 100.0 518.67 642.46 1586.94 1401.34 14.62 ... 2388.03 8130.10 8.4441 0.03 393 2388 100.0 39.08 23.4166 140
3 1 4 0.0042 0.0000 100.0 518.67 642.44 1584.12 1406.42 14.62 ... 2388.05 8132.90 8.3917 0.03 391 2388 100.0 39.00 23.3737 139
4 1 5 0.0014 0.0000 100.0 518.67 642.51 1587.19 1401.92 14.62 ... 2388.03 8129.54 8.4031 0.03 390 2388 100.0 38.99 23.4130 138

5 rows × 27 columns

4.3.2. Normalize¶

In [18]:
norm_test_df = pd.DataFrame(min_max_scaler.transform(test_df[cols_normalize]), 
                            columns=cols_normalize, 
                            index=test_df.index)

test_join_df = test_df[test_df.columns.difference(cols_normalize)].join(norm_test_df)
test_df = test_join_df.reindex(columns = test_df.columns)
test_df = test_df.reset_index(drop=True)

test_df.head()
Out[18]:
id cycle setting1 setting2 setting3 s1 s2 s3 s4 s5 ... s13 s14 s15 s16 s17 s18 s19 s20 s21 RUL
0 1 1 0.632184 0.750000 0.0 0.0 0.545181 0.310661 0.269413 0.0 ... 0.220588 0.132160 0.308965 0.0 0.333333 0.0 0.0 0.558140 0.661834 142
1 1 2 0.344828 0.250000 0.0 0.0 0.150602 0.379551 0.222316 0.0 ... 0.264706 0.204768 0.213159 0.0 0.416667 0.0 0.0 0.682171 0.686827 141
2 1 3 0.517241 0.583333 0.0 0.0 0.376506 0.346632 0.322248 0.0 ... 0.220588 0.155640 0.458638 0.0 0.416667 0.0 0.0 0.728682 0.721348 140
3 1 4 0.741379 0.500000 0.0 0.0 0.370482 0.285154 0.408001 0.0 ... 0.250000 0.170090 0.257022 0.0 0.250000 0.0 0.0 0.666667 0.662110 139
4 1 5 0.580460 0.500000 0.0 0.0 0.391566 0.352082 0.332039 0.0 ... 0.220588 0.152751 0.300885 0.0 0.166667 0.0 0.0 0.658915 0.716377 138

5 rows × 27 columns

4.4. Data Generation¶

In [19]:
# function to generate input sequences
def gen_sequence(id_df, seq_length, seq_cols):
    data_matrix = id_df[seq_cols].values
    num_elements = data_matrix.shape[0]
    
    for start, stop in zip(range(0, num_elements-seq_length), range(seq_length, num_elements)):
        yield data_matrix[start:stop, :]
        
# function to generate labels
def gen_labels(id_df, seq_length, label):
    data_matrix = id_df[label].values
    num_elements = data_matrix.shape[0]
    
    return data_matrix[seq_length:num_elements, :]
In [20]:
# generate data
n_step = 50

sensor_cols = ['s' + str(i) for i in range(1,22)]
x_cols = ['setting1', 'setting2', 'setting3']
x_cols.extend(sensor_cols)

X_train = np.concatenate(list(list(gen_sequence(train_df[train_df['id']==id], n_step, x_cols)) for id in train_df['id'].unique()))
y_train = np.concatenate(list(list(gen_labels(train_df[train_df['id']==id], n_step, ['RUL']) for id in train_df['id'].unique())))

print(X_train.shape)
print(y_train.shape)
(15631, 50, 24)
(15631, 1)

5. Build a Model¶

In [21]:
n_step = 50
n_input = train_df.shape[1] - 3

## LSTM shape
n_lstm1 = 100
n_lstm2 = 50

## Fully connected
n_hidden = 100
n_output = 1

n_batch = 64
In [22]:
weights = {
    'hidden' : tf.Variable(tf.random_normal([n_lstm2, n_hidden], stddev = 0.01)),
    'output' : tf.Variable(tf.random_normal([n_hidden, n_output], stddev = 0.01))
}

biases = {
    'hidden' : tf.Variable(tf.random_normal([n_hidden], stddev = 0.01)),
    'output' : tf.Variable(tf.random_normal([n_output], stddev = 0.01))
}

x = tf.placeholder(tf.float32, [None, n_step, n_input])
y = tf.placeholder(tf.float32, [None, n_output])
In [23]:
def LSTM_model(x, weights, biases):
    with tf.variable_scope('rnn'):

        with tf.variable_scope('lstm1'):            
            lstm1 = tf.nn.rnn_cell.LSTMCell(n_lstm1)
            h1, c1 = tf.nn.dynamic_rnn(lstm1, x, dtype = tf.float32)
        with tf.variable_scope('lstm2'):            
            lstm2 = tf.nn.rnn_cell.LSTMCell(n_lstm2)
            h2, c2 = tf.nn.dynamic_rnn(lstm2, h1, dtype = tf.float32)
            
        # Build classifier
        hidden = tf.add(tf.matmul(h2[:,-1,:], weights['hidden']), biases['hidden'])
        hidden = tf.nn.relu(hidden)
        output = tf.add(tf.matmul(hidden, weights['output']), biases['output'])
        return output
In [24]:
LR = 0.001

pred = LSTM_model(x, weights, biases)
loss = tf.square(tf.subtract(y, pred))
loss = tf.reduce_mean(loss)

optm = tf.train.AdamOptimizer(LR).minimize(loss)
WARNING: Logging before flag parsing goes to stderr.
W0118 18:06:11.669138 34360 deprecation.py:323] From <ipython-input-23-3a577a40a0e4>:5: LSTMCell.__init__ (from tensorflow.python.ops.rnn_cell_impl) is deprecated and will be removed in a future version.
Instructions for updating:
This class is equivalent as tf.keras.layers.LSTMCell, and will be replaced by that in Tensorflow 2.0.
W0118 18:06:11.671089 34360 deprecation.py:323] From <ipython-input-23-3a577a40a0e4>:6: dynamic_rnn (from tensorflow.python.ops.rnn) is deprecated and will be removed in a future version.
Instructions for updating:
Please use `keras.layers.RNN(cell)`, which is equivalent to this API
W0118 18:06:11.722839 34360 deprecation.py:506] From c:\users\juwonna7\appdata\local\programs\python\python36\lib\site-packages\tensorflow\python\ops\init_ops.py:1251: calling VarianceScaling.__init__ (from tensorflow.python.ops.init_ops) with dtype is deprecated and will be removed in a future version.
Instructions for updating:
Call initializer instance with the dtype argument instead of passing it to the constructor
W0118 18:06:11.729649 34360 deprecation.py:506] From c:\users\juwonna7\appdata\local\programs\python\python36\lib\site-packages\tensorflow\python\ops\rnn_cell_impl.py:961: calling Zeros.__init__ (from tensorflow.python.ops.init_ops) with dtype is deprecated and will be removed in a future version.
Instructions for updating:
Call initializer instance with the dtype argument instead of passing it to the constructor

6. Batch Maker¶

In [25]:
dataset_train = tf.data.Dataset.from_tensor_slices((X_train, y_train))
dataset_train = dataset_train.shuffle(100000).repeat().batch(n_batch)
iterator_train = dataset_train.make_one_shot_iterator()
next_batch_train = iterator_train.get_next()
W0118 18:06:13.782006 34360 deprecation.py:323] From <ipython-input-25-cbb778ae917e>:3: DatasetV1.make_one_shot_iterator (from tensorflow.python.data.ops.dataset_ops) is deprecated and will be removed in a future version.
Instructions for updating:
Use `for ... in dataset:` to iterate over a dataset. If using `tf.estimator`, return the `Dataset` object directly from your input function. As a last resort, you can use `tf.compat.v1.data.make_one_shot_iterator(dataset)`.

7. Optimize¶

In [26]:
sess = tf.Session()
init = tf.global_variables_initializer()
sess.run(init)

n_iter = 5000
n_prt = 250

for i in range(n_iter):
    train_x, train_y = sess.run(next_batch_train)
    
    sess.run(optm, feed_dict = {x: train_x, y: train_y.reshape(-1,1)})
    c = sess.run(loss, feed_dict = {x: train_x, y: train_y.reshape(-1,1)})
    
    if i % n_prt == 0:
        print("Iter : {}".format(i), "Cost : {}".format(c))
Iter : 0 Cost : 10177.833984375
Iter : 250 Cost : 2810.477783203125
Iter : 500 Cost : 3452.3369140625
Iter : 750 Cost : 1818.179443359375
Iter : 1000 Cost : 640.7874755859375
Iter : 1250 Cost : 862.000732421875
Iter : 1500 Cost : 380.43408203125
Iter : 1750 Cost : 574.8950805664062
Iter : 2000 Cost : 501.2654724121094
Iter : 2250 Cost : 662.8279418945312
Iter : 2500 Cost : 412.4056701660156
Iter : 2750 Cost : 567.5744018554688
Iter : 3000 Cost : 634.7484130859375
Iter : 3250 Cost : 340.80670166015625
Iter : 3500 Cost : 546.6911010742188
Iter : 3750 Cost : 615.843017578125
Iter : 4000 Cost : 384.20947265625
Iter : 4250 Cost : 260.69677734375
Iter : 4500 Cost : 395.52911376953125
Iter : 4750 Cost : 510.26361083984375

8. Test or Evaluation¶

In [27]:
test_x = [test_df[test_df['id']==id][x_cols].values[-n_step:]
                       for id in test_df['id'].unique() if len(test_df[test_df['id']==id]) >= n_step]

test_x = np.asarray(test_x).astype(np.float32)
print(test_x.shape)
(93, 50, 24)
In [28]:
y_mask = [len(test_df[test_df['id']==id]) >= n_step for id in test_df['id'].unique()] 

test_y = test_df.groupby('id')['RUL'].nth(-1)[y_mask].values
test_y = test_y.reshape(test_y.shape[0],1).astype(np.float32)
print(test_y.shape)
(93, 1)
In [29]:
test_pred = sess.run(pred, feed_dict = {x: test_x})

plt.figure(figsize=(15,8))
plt.plot(test_y.ravel(), label = 'Real')
plt.plot(test_pred.ravel(), label = 'Prediction')
plt.legend(fontsize = 15)
plt.show()